/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/operation.hxx"

#include "linal/device/host/operation.hxx"

namespace linal {

  static int apply_pivots_ln  ( Hier_ piv, Hier_ A );
  static int apply_pivots_lt  ( Hier_ piv, Hier_ A );
  static int apply_pivots_rn  ( Hier_ piv, Hier_ A );
  static int apply_pivots_rt  ( Hier_ piv, Hier_ A );

  /*!
    Apply permutation of p to matrix A : left + no trans
  */
  int apply_pivots_ln(Hier_ p, Hier_ A) {
    FLA_Obj AT, AB, A0, A1, A2;
    FLA_Obj pt, pb, p0, p1, p2;

    FLA_Part_2x1( ~A, &AT,
                      &AB, 0, FLA_TOP );
    FLA_Part_2x1( ~p, &pt,
                      &pb, 0, FLA_TOP );

    while ( FLA_Obj_length( pt ) < FLA_Obj_length( ~p ) ) {
      FLA_Repart_2x1_to_3x1( AT, &A0,
                                 &A1,
                             AB, &A2, 1, FLA_BOTTOM );
      FLA_Repart_2x1_to_3x1( pt, &p0,
                                 &p1,
                             pb, &p2, 1, FLA_BOTTOM );
      //------------------------------------------------------------
      FLA_Obj AA;
      FLA_Merge_2x1( A1,
                     A2, &AA );

      Hier_ HH(AA), p_1(p1);
      for (int i=0;i<HH.get_n();++i) {
        Hier_ T;
        HH.extract( T, HH.get_m(), 1, 0, i );

#pragma omp task firstprivate( p_1, T )
        {
          host::internal::apply_pivots_ln( ((FLA_Obj*)p_1.get_buffer())[0],
                                           ~(T.flat()) );
        }

      }

#pragma omp taskwait

      //------------------------------------------------------------
      FLA_Cont_with_3x1_to_2x1( &AT, A0,
                                A1,
                                &AB, A2, FLA_TOP );
      FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                p1,
                                &pb, p2, FLA_TOP );
    }
    return true;
  }

  /*!
    Apply permutation of p to matrix A : left + trans 
  */
  int apply_pivots_lt(Hier_ p, Hier_ A) {
    FLA_Obj AT, AB, A0, A1, A2;
    FLA_Obj pt, pb, p0, p1, p2;

    FLA_Part_2x1( ~A, &AT,
                  &AB, 0, FLA_BOTTOM );
    FLA_Part_2x1( ~p, &pt,
                  &pb, 0, FLA_BOTTOM );

    while ( FLA_Obj_length( pb ) < FLA_Obj_length( ~p ) ) {
      FLA_Repart_2x1_to_3x1( AT, &A0,
                             &A1,
                             AB, &A2, 1, FLA_TOP );
      FLA_Repart_2x1_to_3x1( pt, &p0,
                             &p1,
                             pb, &p2, 1, FLA_TOP );
      //------------------------------------------------------------
      FLA_Obj AA;
      FLA_Merge_2x1( A1,
                     A2, &AA );

      Hier_ HH(AA), p_1(p1);
      for (int i=0;i<HH.get_n();++i) {
        Hier_ T;
        HH.extract( T, HH.get_m(), 1, 0, i );

#pragma omp task firstprivate( p_1, T )
        {
          host::internal::apply_pivots_lt( ((FLA_Obj*)p_1.get_buffer())[0],
                                           ~(T.flat()) );
        }

      }

#pragma omp taskwait

      //------------------------------------------------------------
      FLA_Cont_with_3x1_to_2x1( &AT, A0,
                                A1,
                                &AB, A2, FLA_BOTTOM );
      FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                p1,
                                &pb, p2, FLA_BOTTOM );
    }
    return true;
  }


  int apply_pivots_rn(Hier_ p, Hier_ A) {
    FLA_Obj AL, AR, A0, A1, A2;
    FLA_Obj pt, pb, p0, p1, p2;

    FLA_Part_1x2( ~A, &AL, &AR, 0, FLA_RIGHT );
    FLA_Part_2x1( ~p, &pt,
                  &pb, 0, FLA_BOTTOM );

    while ( FLA_Obj_length( pb ) < FLA_Obj_length( ~p ) ) {
      FLA_Repart_1x2_to_1x3( AL, AR, &A0, &A1, &A2, 1, FLA_LEFT );
      FLA_Repart_2x1_to_3x1( pt, &p0,
                             &p1,
                             pb, &p2, 1, FLA_TOP );
      //------------------------------------------------------------
      FLA_Obj AA;
      FLA_Merge_1x2( A1, A2, &AA );

      Hier_ HH(AA), p_1(p1);
      for (int i=0;i<HH.get_m();++i) {
        Hier_ T;
        HH.extract( T, 1, HH.get_n(),  i, 0 );

#pragma omp task firstprivate( p_1, T )
        {
          host::internal::apply_pivots_rn( ((FLA_Obj*)p_1.get_buffer())[0],
                                           ~(T.flat()) );
        }
      }

#pragma omp taskwait

      //------------------------------------------------------------
      FLA_Cont_with_1x3_to_1x2( &AL, &AR, A0, A1, A2, FLA_RIGHT );
      FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                p1,
                                &pb, p2, FLA_BOTTOM );
    }
    return true;
  }

  int apply_pivots_rt(Hier_ p, Hier_ A) {
    FLA_Obj AL, AR, A0, A1, A2;
    FLA_Obj pt, pb, p0, p1, p2;

    FLA_Part_1x2( ~A, &AL, &AR, 0, FLA_LEFT );
    FLA_Part_2x1( ~p, &pt,
                  &pb, 0, FLA_TOP );

    while ( FLA_Obj_length( pt ) < FLA_Obj_length( ~p ) ) {
      FLA_Repart_1x2_to_1x3( AL, AR, &A0, &A1, &A2, 1, FLA_RIGHT );
      FLA_Repart_2x1_to_3x1( pt, &p0,
                             &p1,
                             pb, &p2, 1, FLA_BOTTOM );
      //------------------------------------------------------------
      FLA_Obj AA;
      FLA_Merge_1x2( A1, A2, &AA );

      Hier_ HH(AA), p_1(p1);
      for (int i=0;i<HH.get_m();++i) {
        Hier_ T;
        HH.extract( T, 1, HH.get_n(),  i, 0 );

#pragma omp task firstprivate( p_1, T )
        {
          host::internal::apply_pivots_rt( ((FLA_Obj*)p_1.get_buffer())[0],
                                           ~(T.flat()) );
        }
      }

#pragma omp taskwait

      //------------------------------------------------------------
      FLA_Cont_with_1x3_to_1x2( &AL, &AR, A0, A1, A2, FLA_LEFT );
      FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                p1,
                                &pb, p2, FLA_TOP );
    }
    return true;
  }

  int apply_pivots(int side, int trans, Hier_ p, Hier_ A) {
    LINAL_ERROR( p.get_datatype() == FLA_INT,
                 LINAL_ERROR_CHECK_INT );
      
    switch (side) {
    case FLA_LEFT: 
      {
        switch (trans) {
        case FLA_NO_TRANSPOSE: 
          apply_pivots_ln(p, A);    
          break;
        case FLA_TRANSPOSE:
          apply_pivots_lt(p, A);
          break;
        }
        break;
      }
    case FLA_RIGHT:
      {
        switch (trans) {
        case FLA_NO_TRANSPOSE:
          apply_pivots_rn(p, A);
          break;
        case FLA_TRANSPOSE:
          apply_pivots_rt(p, A);
          break;
        }
        break;
      }
    default:
      {
        LINAL_ERROR_DETAIL( check_side     ( side  ),
                            "%s\n %s %d\n",
                            LINAL_ERROR_CHECK_SIDE,
                            "side = ", side );
        LINAL_ERROR_DETAIL( check_trans_only( trans ),
                            "%s\n %s %d\n",
                            LINAL_ERROR_CHECK_TRANS,
                            "trans = ", trans );
        break;
      }
    }
    return true;      
  }
}

